rm(list=ls())

#---------------------------------------------------------------#
# PROJECT: HUMAN RIGHTS ABUSES                                  #
# SCRIPT: REPLICATION MANUSCRIPT - FIG 2,3,4,7,8,9              #
# AUTHOR: FLORES MACIAS AND ZARKIN                              #
# DATE: DECEMBER 2022                                           #
#---------------------------------------------------------------#

#### ------ STEP 1: LOAD PACKAGES #### 

library(foreign)
library(tidyverse)
library(readstata13)
library(ggplot2)
library(PanelMatch)

#### ------ STEP 2: DEFINE THEME FOR FIGURES #### 


theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90),
      axis.title.x =      element_text(size =12),
      legend.position = "bottom",
      panel.background = element_rect(fill = "white"),
      strip.background = element_rect(colour="black", fill="white")
    )
}

#### ------ STEP 3: LOAD DATA SETS #### 

hr<- read.dta13("0_DATA/ForMatching.dta")
hra<- read.dta13("0_DATA/ForMatching_Always.dta")
hrd<- read.dta13("0_DATA/ForMatching_NatDisOnly.dta")
  hrd <- subset(hrd, year >= 2010)

#### ------ STEP 4: I, T IN INTEGER FORM #### 

hr$year <- as.integer(hr$year)
hr$inegi <- as.integer(hr$inegi)

hra$year <- as.integer(hra$year)
hra$inegi <- as.integer(hra$inegi)

hrd$year <- as.integer(hrd$year)
hrd$inegi <- as.integer(hrd$inegi)



#### ------ STEP 5: FIGURE 2. UNITS BY TREATMENT AND CONTROL #### 

trfigop<- DisplayTreatment(unit.id = "inegi",
                           time.id = "year", legend.position = "bottom",
                           xlab = "Year", ylab = "Municipality",
                           treatment = "operativo", data = hr,
                           color.of.treated = "#636363",
                           color.of.untreated = "#f0f0f0",
                           legend.labels = c("Not Treated", "Treated"),
                           hide.y.tick.label = TRUE,
                           dense.plot = TRUE,
                           x.angle=90)
trfigop<-trfigop + labs(fill="", title="")
plot(trfigop)

ggsave("2_FIGURES/FIG2.png", plot=trfigop, width=10,height=8, units = "in", dpi = 300)

#### ------ STEP 6: FIGURE 3. COVARIATE BALANCE BY REFINEMENT METHOD #### 

##------------ NO MATCHING METHOD ####

# MATCHING
none_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                         treatment = "operativo", refinement.method = "none", 
                         data = hr, match.missing = T, 
                         covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3))  + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                         qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                         exact.match.variables = c("gradorezagosoc2010num", "rural"),
                         lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(none_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_none<-data.frame(get_covariate_balance(none_match$att, 
                                           hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                           plot = F, ylim = c(-2,2)))

##------------ PS MATCH 5 ####


# MATCHING
psm5_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                         treatment = "operativo", refinement.method = "ps.match", 
                         data = hr, match.missing = T, 
                         covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                         size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                         exact.match.variables = c("gradorezagosoc2010num", "rural"),
                         lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(psm5_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_psm5<-data.frame(get_covariate_balance(psm5_match$att, 
                                           hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                           plot = F, ylim = c(-2,2)))

##------------ PS MATCH 10 ####


#MATCHING
psm10_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.match", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3))  + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                          size.match = 10, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:3, forbid.treatment.reversal = FALSE)

#GET COVARIATE BALANCE
get_covariate_balance(psm10_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_psm10<-data.frame(get_covariate_balance(psm10_match$att, 
                                            hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                            plot = F, ylim = c(-2,2)))

##------------ CBPS MATCH 5 ####


# MATCHING
cbpsm5_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "CBPS.match", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                           size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(cbpsm5_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_cbps5<-data.frame(get_covariate_balance(cbpsm5_match$att, 
                                            hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                            plot = F, ylim = c(-2,2)))
##------------ CBPS MATCH 10 ####


cbpsm10_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.match", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                            size.match = 10, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:3, forbid.treatment.reversal = FALSE)

get_covariate_balance(cbpsm10_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_cbps10<-data.frame(get_covariate_balance(cbpsm10_match$att, 
                                             hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                             plot = F, ylim = c(-2,2)))

##------------ MAHALANOBIS ####


# MATCHING
mahal_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "mahalanobis", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(mahal_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_mahal<-data.frame(get_covariate_balance(mahal_match$att, 
                                            hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                            plot = F, ylim = c(-2,2)))

##------------ PS WEIGHTS ####


# MATCHING
psw_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                        treatment = "operativo", refinement.method = "ps.weight", 
                        data = hr, match.missing = T, 
                        covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                        qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                        exact.match.variables = c("gradorezagosoc2010num", "rural"),
                        lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(psw_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_psw<-data.frame(get_covariate_balance(psw_match$att, 
                                          hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                          plot = F, ylim = c(-2,2)))

##------------ CBPS WEIGHTS ####


cbpsw_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "CBPS.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:3, forbid.treatment.reversal = FALSE)

get_covariate_balance(cbpsw_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
cov_cbpsw<-data.frame(get_covariate_balance(cbpsw_match$att, 
                                            hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samegovpres"), 
                                            plot = F, ylim = c(-2,2)))

##------------ FIGURE ####


# GENERATE NEW VARIABLES AND RENAME

cov_cbps10$Time<-c("1","2","3")
cov_cbps5$Time<-c("1","2","3")
cov_cbpsw$Time<-c("1","2","3")
cov_mahal$Time<-c("1","2","3")
cov_none$Time<-c("1","2","3")
cov_psm10$Time<-c("1","2","3")
cov_psm5$Time<-c("1","2","3")
cov_psw$Time<-c("1","2","3")

cov_cbps10<-cov_cbps10 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_cbps5<-cov_cbps5 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_cbpsw<-cov_cbpsw %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_mahal<-cov_mahal %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_none<-cov_none %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_psm10<-cov_psm10 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_psm5<-cov_psm5 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_psw<-cov_psw %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")


cov_cbps10$Method<-c("CBPS match 10")
cov_cbps5$Method<-c("CBPS match 5")
cov_cbpsw$Method<-c("CBPS weights")
cov_mahal$Method<-c("Mahalanobis")
cov_none$Method<-c("None")
cov_psm10$Method<-c("PS match 10")
cov_psm5$Method<-c("PS match 5")
cov_psw$Method<-c("PS weights")


cov_cbps10$Covariate<-recode(cov_cbps10$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_cbps5$Covariate<-recode(cov_cbps5$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_cbpsw$Covariate<-recode(cov_cbpsw$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_mahal$Covariate<-recode(cov_mahal$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_none$Covariate<-recode(cov_none$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_psm10$Covariate<-recode(cov_psm10$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_psm5$Covariate<-recode(cov_psm5$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")
cov_psw$Covariate<-recode(cov_psw$Covariate, "thom" = "Homicide rate", "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population", "turfwar_threeaf" = "Turf war", "samegovpres" = "Same party")


# APPEND DATA SETS

cov_ops<-rbind(cov_cbps10,cov_cbps5,cov_cbpsw,cov_mahal,
               cov_none, cov_psm10, cov_psm5, cov_psw)

cov_ops$Method <- factor(cov_ops$Method, levels = c("None", "PS weights", "PS match 5", "PS match 10", 
                                                    "CBPS weights", "CBPS match 5", "CBPS match 10", "Mahalanobis"))
cov_ops$Time<-as.numeric(cov_ops$Time)

# CREATE FIGURE

covopsfig=ggplot(cov_ops, aes(x=Time, y=Estimate, color=Covariate)) +
  geom_line(size=1) +
  facet_wrap(vars(Method), nrow=2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="black",  linetype=2) +
  scale_x_continuous(breaks=c(1,2,3), labels=c("t-3", "t-2", "t-1")) +
  scale_y_continuous(name="Mean difference (st. dev)", limits=c(-2,2)) +
  scale_color_grey()

plot(covopsfig)

ggsave("2_FIGURES/FIG3.png", plot=covopsfig, width=12,height=6, units = "in", dpi = 300)

rm(cbpsm10_match,cbpsm5_match,cbpsw_match,cov_cbps10,cov_cbps5,
   cov_cbpsw,cov_mahal,cov_none,cov_psm10,cov_psm5,cov_psw,
   mahal_match,psm10_match,psm5_match,psw_match,none_match)

#### ------ STEP 7: FIGURE 4. ESTIMATED AVERAGE EFFECTS OF CONSTABULARIZATION ON THE RATE OF SERIOUS HR ABUSE COMPLAINTS AGAINST SECURITY FORCES ####

##------------ PS WEIGHTS ####

# MATCHING  

mainps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:3, forbid.treatment.reversal = FALSE)
# DIFF & DIFF

mainps_estim <- PanelEstimate(sets = mainps_match, data = hr)

summary(mainps_estim)

main_ps<-as.data.frame(summary(mainps_estim))

##------------ CBPS WEIGHTS ####


# MATCHING 

maincbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)
# DIFF & DIFF  

maincbps_estim <- PanelEstimate(sets = maincbps_match, data = hr)

summary(maincbps_estim)

main_cbps<-as.data.frame(summary(maincbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

main_ps$Method<-c("PS weights")
main_cbps$Method<-c("CBPS weights")

main_ps<-main_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

main_cbps<-main_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

main_ps$Time<-c("t+0","t+1","t+2","t+3" )
main_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

main<-rbind(main_cbps,main_ps)

# CREATE FIGURE

main_fig=ggplot(main, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(main_fig)

ggsave("2_FIGURES/FIG4.png", plot=main_fig, width=7,height=5, units = "in", dpi = 300)

rm(main_cbps,main_ps,maincbps_estim,maincbps_match,mainps_estim,mainps_match)

#### ------ STEP 8: FIGURE 7. ESTIMATED AVERAGE EFFECTS OF MILITARY DEPLOYMENTS FOR DISASTER RELIEF ON THE RATE OF SERIOUS HR ABUSE COMPLAINTS AGAINST SECURITY FORCES ####

##------------ PS WEIGHTS ####

# MATCHING

natdisps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                             treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                             data = hrd, match.missing = T, 
                             covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(natdisps_match$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

natdisps_results <- PanelEstimate(sets = natdisps_match, data = hrd)

summary(natdisps_results)

nd_ps<-as.data.frame(summary(natdisps_results))

##------------ CBPS WEIGHTS ####

# MATCHING

natdiscbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                               treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                               data = hrd, match.missing = T, 
                               covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)), 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(natdiscbps_match$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

natdiscbps_results <- PanelEstimate(sets = natdiscbps_match, data = hrd)

summary(natdiscbps_results)

nd_cbps<-as.data.frame(summary(natdiscbps_results))


##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

nd_ps$Method<-c("PS weights")
nd_cbps$Method<-c("CBPS weights")

nd_ps<-nd_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

nd_cbps<-nd_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

nd_ps$Time<-c("t+0","t+1","t+2" )
nd_cbps$Time<-c("t+0","t+1","t+2")

# APPEND DATA SETS

natdis<-rbind(nd_cbps,nd_ps)

# CREATE FIGURE

nd_fig=ggplot(natdis, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(nd_fig)

ggsave("2_FIGURES/FIG7.png", plot=nd_fig, width=7,height=5, units = "in", dpi = 300)

rm(natdiscbps_match, natdiscbps_results, natdisps_match, natdisps_results)

#### ------ STEP 9: FIGURE 8. ESTIMATED AVERAGE EFFECTS OF CONSTABULARIZATION ON THE RATE OF HR ABUSE COMPLAINTS FILED AGAINST THE IMSS ####

##------------ PS WEIGHTS ####

# MATCHING  

imssps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:3))  + I(lag(quejasIMSSpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                           qoi = "att" ,outcome.var = "quejasIMSSpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(imssps_match$att, hr, covariates = c("thom", "quejasIMSSpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

imssps_estim <- PanelEstimate(sets = imssps_match, data = hr)

summary(imssps_estim)

imss_ps<-as.data.frame(summary(imssps_estim))


##------------ CBPS WEIGHTS ####


# MATCHING 

imsscbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(quejasIMSSpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "quejasIMSSpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(imsscbps_match$att, hr, covariates = c("thom", "quejasIMSSpc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF  

imsscbps_estim <- PanelEstimate(sets = imsscbps_match, data = hr)

summary(imsscbps_estim)

imss_cbps<-as.data.frame(summary(imsscbps_estim))


##------------ FIGURE ####


# GENERATE NEW VARIABLES AND RENAME

imss_ps$Method<-c("PS weights")
imss_cbps$Method<-c("CBPS weights")

imss_ps<-imss_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

imss_cbps<-imss_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

imss_ps$Time<-c("t+0","t+1","t+2","t+3" )
imss_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

imss<-rbind(imss_cbps,imss_ps)

# CREATE FIGURE

imss_fig=ggplot(imss, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-0.3, 0.3) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(imss_fig)

ggsave("2_FIGURES/FIG8.png", plot=imss_fig, width=7,height=5, units = "in", dpi = 300)

rm(imsscbps_estim, imsscbps_match, imssps_estim, imssps_match)

#### ------ STEP 10: FIGURE 9. ESTIMATED AVERAGE EFFECTS OF CONSTABULARIZATION OVER TIME IN MUNICIPALITIES WITH UNINTERRUPTED MILITARY PRESENCE ####

##------------ PS WEIGHTS ####

# MATCHING

alwaysps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "ps.weight", 
                             data = hra, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3))  + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)) , 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:8, forbid.treatment.reversal = FALSE)

# DIFF & DIFF  

alwaysps_estim <- PanelEstimate(sets = alwaysps_match, data = hra)

summary(alwaysps_estim)

always_ps<-as.data.frame(summary(alwaysps_estim))


##------------ CBPS WEIGHTS ####


# MATCHING

alwayscbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                               treatment = "operativo", refinement.method = "CBPS.weight", 
                               data = hra, match.missing = T, 
                               covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:8, forbid.treatment.reversal = FALSE)

# DIFF & DIFF 

alwayscbps_estim <- PanelEstimate(sets = alwayscbps_match, data = hra)

summary(alwayscbps_estim)

always_cbps<-as.data.frame(summary(alwayscbps_estim))


##------------ FIGURE ####


# GENERATE NEW VARIABLES AND RENAME

always_ps$Method<-c("PS weights")
always_cbps$Method<-c("CBPS weights")

always_ps<-always_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

always_cbps<-always_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

always_ps$Time<-c("t+0","t+1","t+2","t+3","t+4","t+5","t+6","t+7","t+8")
always_cbps$Time<-c("t+0","t+1","t+2","t+3","t+4","t+5","t+6","t+7","t+8")

# APPEND DATA SETS

always<-rbind(always_cbps,always_ps)

# CREATE FIGURE

always_fig=ggplot(always, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-5, 5) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(always_fig)

ggsave("2_FIGURES/FIG9.png", plot=always_fig, width=8,height=4, units = "in", dpi = 300)

rm(alwayscbps_estim, alwayscbps_match, alwaysps_estim, alwaysps_match)


